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The statements in the DANSKY subroutine 

30 A (J.NMI) = SAVE 
DO *10 J = 1, NN 


should read 

(unchanged) 

(changed) 


90 SUM = SUM + A (NMIP1 ,K)*A(K,1 1) 
I F ( I I .LT.NMI) GO TO 100 
I F ( I I .EQ.NN) GO TO 100 
SUM = SUM + A (NM I PI ,11 + 1) 

100 A (NM I , II) = SUM 


(unchanged) 

(unchanged) 

(unchanged) 

(unchanged) 

(unchanged but three 
duplicate lines removed) 


Page 1*»: The statements should read 

CALL BOLLIN (A , AS , B , C , X , Y , Z1 , Z2 , Z3 , N , NMAX) 

DOUBLE PRECISION: AS, X,Y,Z1,Z2,Z3 must be double precision 

NOTES: (l) A and B may be destroyed in call to DAVISO 


This information is being published in prelimi- 
nary form in order to expedite its early release. 
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ABSTRACT 


FORTRAN computer subroutines stemming from requirements to pro- 
cess state variable system equations for systems of high order are 
presented. They find the characteristic equation of a matrix using 
the method of Danilevsky, the number of roots with positive real 
parts using the Routh-Horwit 2 alternate formulation, convert a state 
variable system description to a Laplace transfer function using 
the method of Bollinger, and evaluate that transfer function and 
obtain its frequency response. A sample problem is presented to 
demonstrate use of the subroutines. 
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COMPUTER PROGRAMS FOR CALCULATION OF MATRIX STABILITY AND 


FREQUENCY RESPONSE FROM A STATE-SPACE SYSTEM DESCRIPTION 

by Robert C. Seidel 
Lewis Research Center 


SUMMARY 

FORTRAN computer subroutines stemming from requirements to pro- 
cess state variable system equations for systems of high order are 
presented. They find the characteristic equation of a matrix using 
the method of Danilevsky, the number of roots with positive real 
parts using the Routh-Horwitz alternate formulation, convert a state 
variable system description to a Laplace transfer function using the 
method of Bollinger, and evaluate that transfer function and obtain 
its frequency response. A sample problem is presented to demonstrate 
use of the subroutines. 


INTRODUCTION 


High-speed digital computers have made matrix state variable 
methods for system analysis practical. But for large systems, the 
required execution time and the cumulative effect of round off errors 
make it increasingly important to employ efficient algorithms. 

FORTRAN subroutines resulting from requirements to handle large 
systems are reported herein. In particular, they find the charac- 
teristic equation of a system matrix, test that equation for the num- 
ber of roots with positive real parts, convert a state variable 
system description to a Laplace transfer function and evaluate the 
transfer function at a given frequency to obtain its frequency response. 

The FORTRAN program for obtaining the characteristic equation of 
a system matrix uses the method of Danilevsky, reference 1. The pro- 
gram includes a Gauss pivotal element condensation scheme to somewhat 
increase its accuracy. A competing method is that of Leverrier, 
references 1 and 2. However, to compute the characteristic polynomial 
the Danilevsky method is known to be more accurate, faster, and re- 
quire less storage. The FORTRAN program for determining the stability 
of the characteristic polynomial uses the Routh-Horwitz alternate 
formulation method, reference 3. The FORTRAN program for obtaining 
the system transfer function uses the method of Bollinger, reference 4. 
This method appears more reliable tl^an the method of Davison, refer- 
ence 5, which in certain cases is known to converge improperly 
(ref. 4). However Davison's transformation, which permits the output 
to be an arbitrary linear combination of the states, is used in 
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conjunction with Bollinger’s method. The inverse of this transfor- 
mation is also required, but is easily calculated in closed form. 
Once the transfer function is obtained its frequency response is 
easily calculated. 

The FORTRAN listings and subroutine descriptions are presented 
in the following section. Equation symbols are defined in Appendix 
A. FORTRAN symbols are defined separately for each program. A 
sample problem demonstrating the use of the subroutines is described 
in Appendix B. 


COMPUTER PROGRAMS 


Characteristic Equation 

In subroutine DANSKY the characteristic equation of a matrix is 
found by the method of Danilevsky, reference 1. The characteristic 
equation of an nxn matrix A is an expansion of the determinant 
equation 


| A - XI | = <-l) n (X n + Zl^ 11 " 1 + ••• + Zl n ) = 0 

where I Is the identity matrix and the polynomial coefficients 
sought are in the Z1 vector. The Z1 vector coefficients are 
obtained in DANSKY through successive application of similarity 
transformations which finally produce the Zl vector in the top row 
of A. As noted in reference 1 the method allows use of a Gauss 
pivotal element scheme to somewhat increase its accuracy. In DANSKY 
this option is implemented. The Gauss method performs similarity 
transformations which interchange columns and rows of A to place the 
element with the largest absolute value in pivot position. The 
execution time for an 18x18 A matrix is about 0.35 sec (IBM-360-67 TSS 
computer) . DANSKY uses (as most of the programs) double precision. 

The time penalty for using double precision is only about 0.02 sec 
for the 18x18 matrix. One problem which occurs with certain A matrices 
when using DANSKY is that of exponent under or overflow. In such 
cases the matrix A may be (time) scaled by multiplying each element 
of A by a positive constant, r. The characteristic polynomial of the 
scaled matrix becomes 

X n + (Zl^JX 11 " 1 + ••• + (Zl n r n ) 

Figure 1 Is a description of the DANSKY calling statement transfer 
variables. Figure 2 is a FORTRAN listing of DANSKY. 
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For comparison with the Danilevsky method some results obtained 
with the Leverrier method are cited. For the same 18x18 matrix dis- 
cussed earlier, the Leverrier program required 4.2 sec (compared to 
0.35 sec for the Danilevsky method). Also, two additional 18x18 
double precision scratch storage matrices are required for the 
Leverrier method, and double precision is required to obtain the same 
accuracy achieved by DANSKY in single precision. 

In DANSKY a call is made to subroutine POLMPY to multiply two 
polynomials. In certain cases the method of Danilevsky obtains the 
characteristic equation coefficients in partially factored form and 
the factors must be multiplied together to obtain the characteristic 
equation. In DANSKY a call is made to POLMPY in all cases with one 
of the polynomials possibly unity. Figure 3 is a description of the 
POLMPY transfer variables. Figure 4 is a FORTRAN listing of POLMPY, 


Stability 

The subroutine RHWTZ performs a stability test upon the charac- 
teristic equation. It counts the number of roots with positive real 
parts without actually finding them, A simple recursive algorithm 
which is well suited to machine computation is used. Its description 
is given in reference 3. The program is specialized in that it assumes 
the leading polynomial coefficient is unity as in the form returned 
by DANSKY. The execution time for an eighteenth order polynomial is 
about 0.007 sec. Figure 5 presents a description of the RHWTZ trans- 
fer variables and figure 6 presents a FORTRAN listing of RHWTZ. 

For certain equations the test may fail if during the algorithm 
execution a zero appears as a divisor term. In this case the number 
of unstable roots is set to -1 and the message "Test Failed M set to 
-1" is output. 


Transfer Function 

BOLLIN is a subroutine for converting a state variable matrix 
differential equation into an equivalent Laplace transfer function. 
The system equations considered are 


x = Ax + Bu, y 


= C 


where x, B, and C are n vectors, A is an nxn matrix, and u and 
y are input and output scalars. The following steps are taken to 
obtain the system transfer function: 
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1. A call to subroutine DAVISO transforms the A and B system 
matrices using the C vector so as to make the output y a state 
variable in the transformed system, reference 5. The transformed 
system is 


A * = TAT” 1 , B* - TB, and x* = Tx 

such that 

A * A 

x = A x + B u 

where T is the ^dentity matrix except that the MC-th row is over- 
written by the C vector. The integer MC is the position of the 
element of C with the maximum absolute value. The T~1 matrix is 
similar to inverses encountered in the proof of the Danilevsky method, 
reference 1, and can be written down explicitly. The output j of 
the modified system is the MC-th modified state variable, 

2. The denominator polynomial Z1 (characteristic equation) of 
t£ie system transfer function is obtained by a call to DANSKY using the 
A matrix. 

3. The numerator polynomial is obtained in two more calls to 
DANSKY using^Bollinger *s methog, reference 4. First, with the MC-th 
column of A replaced by -B to obtain Z2; then, with a matrix^of 
order (n-1) obtained by deleting the MC-th row and column from A 

to obtain Z3. 

4. The numerator polynomial is computed as Z4 = Z2 - sZ3« 

5. The system transfer function y(s)/u(s) is Z4(s)/Zl(s)» 

Execution times found for various order systems were about 1,0 sec 
for an eighteenth order, 2.7 sec for a twenty-sixth order, and 9.7 sec 
for a forty-first order system. Figure 7 is a description of the 
BOLLIN transfer variables. Figure 8 is a FORTRAN listing of BOLLIN. 
Figure 9 is a description of the DAVISO transfer variables and figure 10 
is a FORTRAN listing of DAVISO. If the A matrix is time scaled by 
multiplication by scalar r, the B vector and frequencies used to 
evaluate the transfer function should also be scaled by r« To handle 
the more general case of multiple input, multiple output systems where 
B and C are matrices, the program calling BOLLIN would start the 
B and C matrices at the appropriate columns in the transfer variable 
list to obtain the desired input/output relation, reference 5. 
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Frequency Response 

The subroutine FRPOLY may be used to evaluate the Laplace trans- 
form Z4(s)/Zl(s) ratio of polynomials for a given frequency s = ju>. 
The evaluation is performed in double precision with the real and 
imaginary parts of the powers of jw handled separately. Figure 11 
is a description of the FRPOLY transfer variables and figure 12 is a 
FORTRAN listing of FRPOLY. To evaluate a system frequency response 
over a range of frequencies, multiple calls to FRPOLY would be made. 


CONCLUDING REMARKS 

Six subroutines % DANSKY, POLMPY, RHWTZ, BOLLIN , DAVISO, and 
FRPOLY were presented and discussed. The chart shown in figure 13 
summarizes the flow from a state variable representation of a system 
through the various subroutines. The major Input and output relations 
and alternate paths for various uses are noted. 
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APPENDIX A 


SYMBOLS 


A 

B 

C 

I 

3 


system matrix, nxn 
input vector, n 
output vector, n 
identity matrix, nxn 


imaginary. 


MC position of element in C with maximum absolute value 
n system order 

r scaling factor 

s Laplace variable, sec ^ 

t time, sec 

T transformation matrix, nxn 

u input 

x state variable vector, n 

y output 

Z1 characteristic equation polynomial 

Z2 intermediate polynomial 

Z3 intermediate polynomial 

Z4 system numerator transfer function polynomial 

X root of characteristic equation 

0 ) frequency, hertz 

Superscripts 

denotes transformed variable 
transpose 
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APPENDIX B 


SAMPLE PROBLEM 

A third order sample problem demonstrating the combined use of the 
subroutines is described next. The problem studied is the transfer 
function 


y (s) s + 6 o s + 6 

u(s) (s 2 + 3s + 9)(s + 4) s 3 + 7s 2 + 21s + 36 

In phase variable form the state matrices are 


A - 

0 10 
0 0 1 

, B = 

0 

0 

’ C " 

6 

1 


-36 -21 -7 

L J 


1 


0 


Figure 14 is a FORTRAN listing of the sample problem MAIN program. 
First the A, B, and C matrices are output. Next after a call to B0LLIN, 
the denominator and numerator polynomial coefficients are output. Then, 
after a call to FRPOLY, the transfer function evaluated at one hertz is 
output, and finally, after a call to RHWTZ, the number of unstable roots 
is output. 

Figure 15 is a listing of the program output. The denominator 
polynomial is s J + 7s^ + 21s + 36 with the unity coefficient of s J 
understood. The numerator polynomial is -2.220E-16 s + s + 6. The 
coefficient of the s^ term should actually be zero but due to limited 
numerical precision is slightly in error. The transfer function is 
evaluated at one hertz (s = j2JI) and has a real part of -.03048 and 
an imaginary part of -.01142 or an amplitude of .03255 at -159.5 
degrees. There are no roots of the denominator polynomial with positive 
real parts. 

If only a test of system stability were desired then it would not 
be necessary to use BOLLIN. Instead DANSKY could be called to obtain 
the characteristic equation followed by a call to RHWTZ. BOLLIN is 
organized such that it computes the denominator transfer function Z1 
each time it is called. To save computations the user may wish to modify 
this if the transformed matrix A is known to remain constant for a 
particular set of B and/or C changes. 
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DANSKY 


Danilevsky 

PURPOSE ; To obtain the characteristic equation of a 

square real matrix 


USAGE ; CALL DANSKY (A, X, Y, Z» N, NMAX) 

Where A - square real matrix, order NxN 
X - scratch vector s order N 
Y - scratch vector, order N 

Z - characteristic polynomial vector, order N 
N - order of A 

NMAX - dimension of A in calling program >_ N 

DOUBLE DIMENSION; A must be double dimen- 
sioned (NMAX, NMAX) in 
calling program 

DOUBLE PRECISION; A, X, Y, and Z must be 

double precision in calling 
program 

NOTES; (1) A is destroyed in obtaining Z 

(2) The Z characteristic equation is 
returned in the form 

A n + Z.X 11-1 +■><>»+ Z + Z 
1 n-1 n 

Subroutines called: (1) POLMPY 


FIGURE 1. -DESCRIPTION OF SUBROUTINE DANSKY TRANSFER VARIABLES 
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C COMPUTE THE COEFFICIENTS OF . THF CHARACTER I ST I 0 ECU AT ION 
SUBROUTINE DAMSKYC A, X, Y, Z, N, UMAX > 

DIMENSION A ( UMAX, 1),XC1),Z(1),Y(1> 

DOUBLE PRECISION SUM, A, SAVE, PI VOT,CK, X, Y, Z 
NM=N 

10 NM1-MN-1 

IF(NMl.EQ.O) GO TO 125 
DO 120 1=1, NM1 
NMI-HN-I 
I PVT=NHI 
NMIM1-NMI-1 
NMI P1=NMI +1 

C FIND MAXIMUM ELEMENT IN PIVOT ROW 
XMAX=A(NMI PI, HIM ) 

1 F (NI'IMl. DT. 0) GO TO 15 
I F ( XMAX . Ffl. 0 , ) GO TO HO 
no TO 50 

15 DO 20 Kl=l, NM| Ml 
K=NM I — K1 

I F (ABS( SMDL ( A (NM I PI, K) ) ) , LE • AB8 ( XUAX ) ) 00 20 

XMAX= A ( MM I PI , K ) 

I PVT=K 

20 CONTINUE 

I F ( XMAX . EO. 0 , ) GO TO 140 
I F( I PVT. ED. MM I ) GO TO 50 

C SIMILARITY TRANSFORM SO PIVOT ELEMENT IS THE MAXIMUM 
DO 30 U»1,MMIP1 
SAVE=A( J, I PVT) 

A ( J , I PVT ) =A ( J , NM I ) 

30 A(J,NMI )»SAVE 

DO 1*0 JJ1, MM 
SAVE® A ( I PVT, J) 

A(IPVT,J)=A(NMI,U) 

40 A(NMI,J)»SAVE 

C A UPDATE FOUALS A * MtSUB N-l) 

50 PIV0T«1./A(NM1P1,NMI ) 

DO 80 K»1,NN 
IFCK.EO.NM1) DO TO 80 
CK»-A( NMI PI, K) *PI VOT 
DO 70 1 1 ■!, NMI 

70 Ad l,K)»AO t,K)+A( H,NMI)*RK 

80 CONTINUE 

DO 85 II ® 1 , NM I 

85 A ( 1 1 ,NM1 )«A( 1 1 ,W1I )*PI VDT 

C MULT MtSUR N-l) INVERSE TIMER A UPRATED 
DO 110 1 1 »!, MM 
SUM»fl, 

DO §0 

90 SUM B SUM*A (NMI P1,K)*A(K, IU 

I E ( 1 1 , LT,NMI ) np TO 160 
IFC II ,5D,NN) DO TD 100 
§UM»§UM*A(NMI PI, 1 1 +1 ) 

IFMI.LT.NMI) DD TO 106 
I E ( l ! . 10, NN) DO TD 169 
SUM B §UM+A(NMI PI, 11*1) 

160 A (NMI , II ) b §UM 

110 CONTINUE 

1?0 DQNTINUE 

m NMI P1 B 1 

6 ELEMENTS OF R SIT lOUAL TD WIT ILiMiNTS IN RAN HfM P4 OR =ft 
146 BD HQ d = NMI PI, NN 

NXsU=NMIPl*l 
150 X(NX) S =A(NMI Pl,d) 

0 MULTIPLY DUT PADT0PS DP OUARADTiRI STIR WTIPN 
119 N Y=N«NN 

DD 299 U=M!Y 
260 Y(U)3? 2 (d) 

BALL PDLMPYUlfY,?, NX,NY, NYPNX ) 

9 R1BUPE SYSTEM RRRER 
Nf|sNN-NX 

IFMN.BT.fl) DR TD 19 

RETURN 

iUD 


iiiuss a.^feaTRM lliilsd m luiaeurisi mm 
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POLMPY 


Polynomial Multiplication 

PURPOSE ; To multiply two polynomials X*Y 

CALL POLMPY (X, Y, Z, NX, NY 9 NZ) 

Where X - first polynomial, vector 
Y - second polynomial , vector 
Z - returned product polynomial, vector 
NX - order of X 
NY - order of Y 

NZ = NX + NY (order of Z) returned 

NOTE; (1) leading coefficient of highest 

power Is assumed to be one, l.e,» 
for X; 

.NX . „ -NX— 1 , , 

A + X^A + •• • + 

DOUBLE PRECISION: X, Y, and Z must be 

double precision in the 
calling program 



FIGURE 3. -DESCRIPTION OF SUBROUTINE POLMPY TRANSFER VARIABLES 
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C MULTIPLY X*Y-Z (LEA"!NG POLYNOMIAL COEFF I C I ENTS ASSUMED UNITY) 
SUBROUTINE POI..MPY ( X , Y, Z, NX , NY, NZ ) 

DIMENSION X(l), Y(1),Z(1) 

DOUBLE PRECISION X,Y,Z,YJ 
NZ-NX+MY 

C IF NX = 0 MEANS POI. Y X= 1 jTHERE FORE Z=v 
IF(NX.OT.O) GO TO 6 
DO 4 J-1,NY 
4 Z(J)-Y(J) 

GO TO 40 

C IF NY<=0 MEANS POLY Y=l;THERFEORE Z=X 
6 IF(NY.GT.O) GO TO 15 

DO 8 J = 1 , N X 
8 Z(J)=X(J) 

GO TO 40 

C START MULT I PL I OAT I ON RY MAXING Z = y*s**NX 
15 DO 20 J = 1 , NZ 

YJ = Y( J ) 

IF(J.GT.NY) YJ=0. 

20 Z(U)=YJ 

C MULTI PL I CAT I Oil LOOP 
DO 30 K-1,NX 
Z(K)=Z(K)+X(K) 

DO 25 J=l, NY 
KPJ = K+«I 

25 Z(KPJ)°Z(KP.J)+Y(J)*X(K) 

30 CONTINUE 

40 RETURN 

END 


FIGURE 4. -FORT RAN LISTING FOR SUBROUTINE POLMPY 
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RHWTZ 


ROUTH - HURWITZ ALTERNATE FORMULATION 

To compute the number of roots with positive real 
parts of a polynomial equation 

CALL RHWTZ (C,N,M) 

C - coefficients of polynomial equation, vector 
N - equation order 

M - number of roots with positive real parts 

(M set to -1 if test fails due to attempted 
division by zero) 

The C vector starts with the second coefficient, 
as the first is assumed unity; that is: 

X n + C. A n_1 + • • • + C -X + C =0 
1 n-1 n 

NOTE: C is destroyed upon return 

DOUBLE PRECISION: C must be double precision 

in the calling program 

PROGRAM OUTPUT: "TEST FAILED M SET TO -1" out- 

put if division by zero is 
attempted 

FIGURE 5. -DESCRIPTION OF SUBROUTINE RHWTZ TRANSFER VARIABLES 


PURPOSE : 

USAGE : 

Where 

RESTRICTION: 
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C ROIITH HORUITZ ALTERNATE FORMULATION STABILITY TEST 
SUBROUTINE RH1»TZ<C,N,M) 

DOUBLE PRECISION C, CNM1, COEF, CNM2 
DIMENSION 0(1) 

NMl-N-1 

M-0 

C PROVIDE FOR CASE N=0,1,OR 2 
IF(M.LE.O) 00 TO 30 
I F( N-2 ) it, 8, 10 
4 IF(C(l).LT.O.) M-l 

00 TO 30 
8 CNMl-1. 

00 TO 25 

10 IE(C(l).EQ.O.) GO TO 35 

C0EF-1./C(1) 

C START ALGORITHM LOOP 
no 20 K"2,MMl 
IF(COFF.LT.O) M-M+l 
no 15 J-K, NM1, 2 
15 C(J)=C(J)-C0EF*0( J + l) 

IF(C(K) .EO.O.) GO TO 35 
C0EF=C(K-1)/C(K) 

20 CONTINUE 

C FINISH REMAINING 2 NO ORDER POLYNOMIAL 
CNM2-CCN-2) 

25 I F(CNM2*C( NM1 ) , LT , 0 , ) M-M+l 

23 1 F ( C ( NMD *C( N ) , LT . 0 . ) M-M+l 

30 RETURN 

35 WRITE(6,45) 

45 F0RMAT(1X,23HTEST FAILED M SET TO -1) 

M--1 

GO TO 30 
END 


FIGURE 6. -FORTRAN LISTING FOR SUBROUTINE RHtfTZ 
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Bollinger 

PURPOSE : 

USAGE: 


BOLLIN 


To obtain a system transfer function from system 
matrix equations 

CALL BOLLIN (A,AS,B,X,Y,Z1,Z2,Z3,N,NMAX) 

Where A - system matrix (NxN) 

AS t- scratch matrix (NxN) 

B - system input vector (N) 

C - system output vector (N) 

X - scratch vector (N) 

Y - scratch vector (N) 

Z1 - transfer function denominator, vector (N) 

Z2 transfer function numerator, vector (N) 

Z3 - scratch, vector (N) 

N - system order 

NMAX - dimension of A and AS in calling program N 

DOUBLE DIMENSION: A and AS must be double 

dimensioned (NMAX, NMAX) 
in calling program 

DOUBLE PRECISION: AS must be double precision 

in calling program 

NOTES: (1) A is destroyed in obtaining Z1 and Z2 

(2) The Z1 and Z2 are returned as In 

Z1 = X n + Zl.X n-1 + “•«+Zl 
X n 

Z2 = Z2 X 11 ” 1 + • • »+Z2 

1 n 

SUBROUTINES CALLED : (1) DA VI SO 

(2) DANSKY 


FIGURE 7. -DESCRIPTION OF SUBROUTINE BOLLIN TRANSFER VARIABLES 



C CONVERT X (DOT) ■AX+RU, Y«C ( TRANSPOSE )X TO TRANSFER 
C FUNCTION Y/ll=Z2/Zl RATIO OF POLYNOMIALS 

SUBROUTINE BOLI.I N{ A, AS, B, C, X, Y,Z1, 22 , Z3, N, UMAX) 

01 MENS I ON A(HMAX,l),AS(NMAX,l),B<l),CO.),XU),Y(l> 
OIMENSIOM Z1(1),Z2(1),Z3(1) 

ROUBLE PRECISION AS, X, Y, Zl, Z2, Z3 
C TRANSFORM Z=TX TO MAKP OUTPUT A STATE 
CALL DAVISO(A,B,C,N,NMAX,MC) 

C SAVE A US! NO AS 
DO 20 K«1,N 
00 10 J = 1 , N 
10 AS(K,J)=A(K,J) 

20 CONTINUE 

C FIND REN CHAR EQM COEF’S 

CALL DANSKY(AS,X,Y,Z1,N,NMAX) 

C SET AS=A AGAIN AND OVER WRITE A(K,MC) COLtlfNI WITH -B 
00 40 K-1,N 
00 30 J=1,N 
30 AS(K, J)«A(K,J) 

40 AS(K,MC)=-B(K) 

C FIND FIRST PART OF MUM CHAR EON 
CALL OANSKY ( AS, X, Y, Z 2, N , NMAX ) 

C SET AS=A AGAIN; COLLAPSE MC ROW AND COLUMN 
NM1-N-1 
DO 60 K=1,NM1 
Kl-K 

•IF(K.GE.MC) Kl-K+1 
00 50 J«1,NM1 
Jl-J 

IF(J.GE.MC) Jl-J+1 
50 AS(K,d)»A(Kl,Jl) 

60 CONTINUE 

C FIND SECOND PART OF NUM CHAR EQN 

CALL DANSKY ( AS, X, Y, Z3, NM1, NMAX > 

C SUBSTRACT FIRST PART FROM SECOND PART 
DO 100 J“1,NM1 
100 Z2(J)-Z2(J)-Z3(J) 

RETURN 

END 


FIGURE B.-FORTRAN LISTING FOR SUBROUTINE BOLL IN 
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DAVISO 

Davison 

PURPOSE ; Transform system to make output a state variable 

USAGE ; CALL DAVISO (A,B,C,N,NMAX,MC) 

Where ' A - input system matrix whijh upon return is 
transformed system TAT 

B - input system vector which upon return is 
transformed input vector TB 

C - system output vector 

N - order of A 

NMAX - dimension of A in calling program >_ N 

MC - position of maximum element in C vector 

DOUBLE DIMENSION: A must be double dimensioned 

(NMAX, NMAX) in calling program 

PROGRAM OUTPUT: If all N elements of C are zero 

the message "All elements of C 
are zero" is output and the pro- 
gram is put in PAUSE 

FIGURE 9. -DESCRIPTION OF SUBROUTINE DAVISO TRANSFER VARIABLES 
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C TRANSFORM X(DOT)=AX+BU, Y=C(TRANSPOSE);< USING Z-TY SUCH 
C THAT Y IS A STATE VARIARLE OF Z(DOT )=TAT( INVERSE )+TBU 
SUBROUTINE DAVI SO (A, B, C,N,NMAX, MC) 

D I MENS I OH A ( UMAX, 1 ) , B( 1 ) , C { 1 ) 

DOUBLE PRECISION SUM 
C FIND MAX EI.Et’FMT INCAS C(MC) 

CM-0. 

MC-0 

DO 5 J»1,N 

I F ( AB5(C( J ) ) . LT . CM) 00 TO 5 
CM-ABS(C (J ) ) 

MC-J 

5 CONTINUE 

IF(MC.GT.O) GO TO 10 
WRITE(6,7) 

7 FORMATt IX, 26HALL ELEMENTS OF C ARE ZERO) 

PAUSE 

C PRFMULT A BY T CHANGES MC ROW ONLY 
10 00 .17 K=1,M 

SUM-0. 

DO IS J-l, M 

15 SUM=SUM+DBIECC())*A(J,K)) 

17 A(MC,K)-SUM 

C POST MULT A BY T INVERSE 

C FOR ,1 NOT = MC: A( K, J ) =A( K, J ) -A ( K, MC ) *C ( J ) / C ( MC ) 

PIV0T=1.7C(MC) 

DO 50 J-l , N 
IFCJ.EQ.MC) GO TO 50 
DO 25 K=l, N 

25 A(K,J)=A(K,J)-A(K, MC )*C(J)*PIVOT 

50 CONTINUE 

C FOR J=MC: A{K,MC)-A(K,MC)/C(MC) 

DO 55 K=1,N 

35 A(K,MC)-A(K,MC)*PIYOT 

C PRE MULT B BY T CHANCES BCMC) ONLY 
SUM-0. 

DO 37 J-l, N 

37 SUM-SUM+DBLE (C(J)*B(J) ) 

B(MC)=SUM 

RETURN 

END 


FIGURE 10. -FORTRAN LISTING FOR SUBROUTINE DAVISO 
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Frequency 

PURPOSE : 

USAGE: 


FRPOLY 

Response of Polynomials 

To evaluate system transfer function polynomials 
to obtain system frequency response 

CALL FRPOLY (Z1,Z2,HZ,G,AMP,PHA,N) 

Where Z1 <- denominator polynomial, vector order N 

Z2 t numerator polynomial, vector order N 

HZ - frequency in hertz 

G - system transfer function Z2(jHZ)/Zl(jHZ) 
evaluated at radian frequency HZ*2 tt 

AMP - transfer function amplitude, |g| 

PHA -» transfer function please, i G degrees 

N - ord^r of Zl 

COMPLEX: G must be declared complex in the 

calling program 

DOUBLE PRECISION: Zl apd Z2 must be declared 

double precision in the 
calling program 

NOTES: the Zl and Z2 must be in the form: 

Zl = + Zl.X n,_1 + ••• +Z1 

1 n 

Z2 = Z2- A 11 *" 1 + ... +Z2 

1 n 


FIGURE 11. -DESCRIPTION OF SUBROUTINE FRPOLY TRANSFER VARIABLES 



C EVALUATE TRANSFER FUNCTION Z2(S)/Z1<S) FOP S=G.2P*HZ*.J 
SUBROUTINE FRPOLY(Zl,Z2,HZ,G,AMP, PHA, N) 

DIMENSION Z1(1),Z2(1) 

DOUBLE PREC I S I ON GT,V!J , SUM1R, SUM1 1 , SUM2R, SUM 2 I,ZZ1,ZZ2,Z1,Z2 

COMPLEX C, SUM1, SUM2 

MM1-N-1 

SUMlR-Zl(N) 

SUM1I-0. 

SUM2R“Z2(N) 

SUM2I-0. 

GT-1. 

WJ»HZ*6. 2831853 
1 KK“2 

IF(NMl.EQ.O) GO TO 40 
DO 35 K*l, NM1 
NMK-N-K 
GT-GT*WJ 
ZZ1«Z1(NMK)*GT 
ZZ2-Z2(NMK)*GT 
GO TO <5,10,3,B),KK 
3 ZZ1— ZZ1 

ZZ2--ZZ2 

5 SUM1R-SUM1R+ZZ1 

SUM2R-SDM2R+ZZ2 
GO TO 30 
8 ZZ1--ZZ1 

ZZ2--ZZ2 

10 SUM11-SUM1I+ZZ1 

SUH2 I “SIJM2 I + ZZ2 
30 KK-KK+1 

I FfKK.GT. 4 ) KOI 
35 CONTINUE 

C ADD ON S**N TO DEN SIJM1 

40 SUM1-CMPLX(SN0L(SUM1R),SNGL(SUM1I ) )+GT*WJ*(OMPLX(0.,l.))**(KK- 
SUM2-CMPLX(SNGLCSUM2R),SNGL(SUM2I )) 

G-SUM2/SUM1 
AMP=CABS ( G ) 

PHA-ATAN2(AIMAG{G),REAL(G)) *57. 29578 

RETURN 

END 


FIGURE 12. -FORTRAN LISTING FOR SUBROUTINE FRPOLY j 
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SYSTEM STATE VARIABLE EQUATIONS 

9 > t 

x = Ax + Bu, y = C x 



equation 



frequency 

response 




RHWTZ 

hr" 

stability 


FIGURE 13. -FLOW CHART FOR FORTRAN SUBROUTINES 









C SAMPLE PROBLEM MAIN PROGRAM 

DIMENSION A C 3, 3 ), ASC 3, 3),BC3), C(3),Z1C3),Z2C3),Z3(3),XC3),Y( 3 ) 
DOUBLE PRECISION AS, Zl, Z2, 23, X, Y 
COMPLEX G 

DATA A/0.,0.,-36., 1., 0. , -21. , (1 . , 1 . , -7./ 

DATA B/O , , 0 , , l./,C/6.,l.,0./,HZ,N, NMAX/1. ,3,3/ 

NR I TE ( 6, 1 ) C C AC I , J), J-l, N) , i-l,N) 

1 FORMAT C ' A=',/,IP3E12,3,/,1P3E12.3,/,1P3E12.35 

UR I T E ( 6, 2 ) (B(d),J-l,N) 

2 FORMAT ( ’ B VECTOR- ' , 1P3E12. 3) 

WRITE(G,3) ( C( J) , J=l, N) 

3 FORMAT ( * C VECTOR--', 1P3E12. 3) 

CALL BOLL INCA, AS, B, C, X, Y, Zl, Z2, Z3, N, UMAX) 

WRITEC6, 10) (Z1(J),J«1,N) 

10 FORMAT C 1 DENOMINATOR-', 1P3E12. 3) 

WRITE (6, 20) CZ2(J),d=l,N) 

20 FORMAT ( 1 NUMERATOR- ', 1P3E 12 . 3 ) 

CALL FRPOLY ( Zl, Z 2, HZ, G, AMP, PHA, N) 

WRITEC6, 30) HZ, G, AMP, PHA 

30 FORMAT C ' AT*,F4.1,' HERTZ THE TRANSFER FUNCTION = ',- 
1 1P2E12. 3, /, 1 AMPLITUDE-', 1PE12. 3, ' PHASE- ' , 1PE12 . 3) 

CALL RHWTZC Zl, N, M) 

WR I TEC 6, 4 0 ) M 

40 FORMAT ( ' NUMBER OF UNSTABLE ROOTS-', I 5) 

STOP 

END 


FIGURE 14.- FORTRAN listing for sample problem main program 


A- 



0.000 

1. 000 E 

00 0.000 


0.000 

0.000 

1.000E 00 


-3.600E 01 

-2. 100E 

01 -7.000E 00 

B 

VECTOR- 

0.000 

0.000 1.000E 00 

C 

VECTOR- 

6.000E on 

1.000E 00 0.000 

DENOMINATOR- 

7.000E 

00 2.100E 01 3.600E 


NUMERATOR- -2.220E-16 1.000E 00 6.000E 00 

AT 1.0 HERTZ THE TRANSFER FUNCTION = -3.048E-02 -1.142E-02 

AMPLITUDE- 3. 255E-02 PHASE- -1.S95E 02 

NUMBER OF UNSTABLE ROOTS- 0 


FIGURE 15. -SAMPLE PROBLEM PROGRAM OUTPUT 



